%{
                                                        ARA_Analyze.m V1
                                                     Written by Rajeev Kant
                                                        Sept. 4th, 2020
                                    Coulombe Lab for Cardiovascular Regenerative Engineering
                            Center for Biomedical Engineering, Brown Univeristy, Providence, RI, 02903
                                                
The purpose of this code is to load in an aortic ring assay wholemount
immunofluorescence stain and analyze color channels. It is assumed that nuclear
stains are on the blue channel and Vimentin is on the green channel. To
analyze an image there must be a stacked RGB image as well as individual
images of each stain. See applyRegion.m for instructions on changing input
stains.

The user may need to edit the following lines: 47 (filepath), 54 (image
scale ratio)

This is the primary file to run all analysis from. See specific
subfunctions for their inputs and outputs. The general tasks of this master
function are to:

1. Segment an input RGB histology image into IF and RM regions
2. Apply boundaries to indivdual stains and determine subROIs for analysis
3. Perform nuclear analysis by count, area, nearest neighbors, and orientation within sub-ROIs
4. Perform vimentin network analysis by area in specified sub-ROIs.

Input:
1. Histology RGB file (upload queried at start of run)
2. Scale factor in px/um

Outputs:
1. ARA_QuantSummary.xlsx: A table containing all calculated parameters (see
individual functions for options)
2. ARA
3. All saved images and .mat files from subfunctions.


Original publication for citation: 
Kant, R. J., Bare, C. F., and Coulombe, K. L. K. "Tissues with 
patterned vessels or protein release induce vascular chemotaxis
in an in vitro platform", Tissue Eng. Part A, 2021.
doi:10.1089/ten.TEA.2020.0269

%}

function ARA_Analyze
%% Load in ARA whole mount stain, get scale, and determine analysis boundaries

close all;
[file,path] = uigetfile("C:\Users\USERNAME\Desktop\ARA Analyze\*"); % Change filepath as appropriate
filename = strcat(file(1:end-4));
RGBim = imread(strcat(path,file));

% Prompt user to input px/�m ratio of image, default is 0.8045. Alter as
% needed
scaleflag = 0; % Flag used to determine if upscaling is needed to analysis occurs at the same resolution
strs  = {'0.8045','0.4023'}; % Change as needed
prompt = {'Select image scale in px/�m:'};
scaleAns = listdlg('PromptString',prompt, 'SelectionMode', 'single', 'ListString', strs, 'ListSize',[140,32]);
if scaleAns == 1
    scale = str2double(strs(1,1));
else
    scale = str2double(strs(1,2));
    scaleflag = 1;
end

%% Determine boundaries for cropping and analysis
imshow(RGBim)
disp('Draw circle around aortic ring for analysis and press enter to continue');
ring = drawcircle; % Draw circle to determine ring location
pause;
origin = ring.Center;
radRing = ring.Radius;
close;

RGBguide = insertShape(RGBim,'circle',[origin(:,1), origin(:,2), radRing], 'LineWidth', 10, 'Color','r'); % Ring inner boundary
lines = [origin(:,1), origin(:,2), origin(:,1)+2500/scale, origin(:,2); % Horizontal RM
    origin(:,1), origin(:,2), origin(:,1)-2500/scale, origin(:,2); % Horizontal IF
    origin(:,1), origin(:,2), origin(:,1)+radRing*3*cosd(30), origin(:,2)+radRing*3*sind(30); % RMbot boundary
    origin(:,1), origin(:,2), origin(:,1)+radRing*3*cosd(30), origin(:,2)-radRing*3*sind(30); % RMtop boundary
    origin(:,1), origin(:,2), origin(:,1)+radRing*3*cosd(150), origin(:,2)+radRing*3*sind(150); % IFbot boundary
    origin(:,1), origin(:,2), origin(:,1)+radRing*3*cosd(210), origin(:,2)+radRing*3*sind(210)]; % IFtop boundary

RGBguide = insertShape(RGBguide, 'line', lines, 'LineWidth', 10, 'Color','r');
RGBguide = insertShape(RGBguide,'circle',[origin(:,1), origin(:,2), 3500*scale], 'LineWidth', 10, 'Color','r'); % Ring outer boundary

% Draw polygonal ROIs based on bounding, IF Region should be drawn first, then RM region
disp('Draw IF boundary using guides and doubleclick confirm, then draw RM boundary.');
polyIF = roipoly(RGBguide);
close;
polyR = roipoly(RGBguide);

% Find bounding boxes for cropping regions
boundIF = regionprops(polyIF, 'Area', 'BoundingBox');
polyIF(:,:,2) = polyIF; % Convert ROI mask to RGB image
polyIF(:,:,3) = polyIF(:,:,1);
boundR = regionprops(polyR, 'Area', 'BoundingBox');
polyR(:,:,2) = polyR; % Convert ROI mask to RGB image
polyR(:,:,3) = polyR(:,:,1);
close all;

% Crop initial RGB image based on determined ROI masks
IF = RGBim;
IF(polyIF == 0) = 0;
IF = imcrop(IF,boundIF.BoundingBox);

R = RGBim;
R(polyR == 0) = 0;
R = imcrop(R,boundR.BoundingBox);

if scaleflag == 1  % Scale up region area normalization factor and ring diameter
    radRing = radRing * 2;
    IF = imresize(IF, 2);
    R = imresize(R, 2);
end
% Save files if needed
%save(fullfile(path, 'radRing.mat'),'radRing');
%save(fullfile(path, 'polyIF.mat'),'polyIF');
%save(fullfile(path, 'boundIF.mat'),'boundIF');
%save(fullfile(path, 'polyR.mat'),'polyR');
%save(fullfile(path, 'boundR.mat'),'boundR');
%save(fullfile(path, 'IF.mat'),'IF');
%save(fullfile(path, 'R.mat'),'R');

%% Apply set boundaries to channels of interest - see applyRegion.m

vimFlag = 0;
if isfile(strcat(path, filename,'_Vim.tif')) % Check if Vim file exists based on naming convention
    vimFlag = 1;
else
    disp("No Vimentin stain file detected");
end

if vimFlag == 1  % Run Vim region if applicable
    [GIF, GR, ~] = applyRegion(filename, path, polyIF, boundIF, polyR, boundR, 2, scaleflag, scale);
    %save(fullfile(path, 'GIF.mat'),'GIF');
    %save(fullfile(path, 'GR.mat'),'GR');
else
    disp("No Vimentin stain file detected");
end

[BIF, BR, scale] = applyRegion(filename, path, polyIF, boundIF, polyR, boundR, 3, scaleflag, scale); % Run DAPI channel
%save(fullfile(path, 'BIF.mat'),'BIF');
%save(fullfile(path, 'BR.mat'),'BR');

%% Define subROIs from RGB image of region - see getROIs.m

[IF_pos] = getROIs (IF, scale, 'IF', path);
%save(fullfile(path, 'IF_pos.mat'),'IF_pos');
[R_pos] = getROIs (R, scale, 'R', path);
%save(fullfile(path, 'R_pos.mat'),'R_pos');


%% Segment and analyze nuclei for subROIs within a region - see analyzeNuclei.m
% analyzeNuclei.m additionally has the subfunctions nearestNeighbors.m and
% histProcess.m within

disp('Segmenting and analyzing nuclei...');
[BIF_stats1, BIF_stats2, BIF_stats3, BIF_stats4, BIF_stats5, BIF_stats6, areaVal, nuclei_IF, BIF_sum] = analyzeNuclei(BIF, IF_pos, radRing, scale, 0, 'IF', path, filename);
[BR_stats1, BR_stats2, BR_stats3, BR_stats4, BR_stats5, BR_stats6, ~, nuclei_R, BR_sum] = analyzeNuclei(BR, R_pos, radRing, scale, areaVal, 'R', path, filename);
%save(fullfile(path, 'nuclei_IF.mat'),'nuclei_IF');
%save(fullfile(path, 'nuclei_R.mat'),'nuclei_R');

%% Analyze Vimentin Area for a region - see vimArea.m

if vimFlag == 1 % Run if Vim file exists, otherwise populate summary arrays with 0s for exporting.
    disp('Analyzing vimentin stains...');
    [GIF_sum] = vimArea(GIF, nuclei_IF, IF_pos, path, 'IF');
    [GR_sum] = vimArea(GR, nuclei_R, R_pos, path, 'R');
else
    GIF_sum = zeros(3,10);
    GR_sum = GIF_sum;
end

% Append G channel summary to statistics
BIF_sum = [BIF_sum; GIF_sum];
BR_sum = [BR_sum; GR_sum];

%% Export and save summary to .xlsx

% Save stats to .mat files for processed image
save(fullfile(path, 'BIF_stats1.mat'),'BIF_stats1');
save(fullfile(path, 'BIF_stats2.mat'),'BIF_stats2');
save(fullfile(path, 'BIF_stats3.mat'),'BIF_stats3');
save(fullfile(path, 'BIF_stats4.mat'),'BIF_stats4');
save(fullfile(path, 'BIF_stats5.mat'),'BIF_stats5');
save(fullfile(path, 'BIF_stats6.mat'),'BIF_stats6');

save(fullfile(path, 'BR_stats1.mat'),'BR_stats1');
save(fullfile(path, 'BR_stats2.mat'),'BR_stats2');
save(fullfile(path, 'BR_stats3.mat'),'BR_stats3');
save(fullfile(path, 'BR_stats4.mat'),'BR_stats4');
save(fullfile(path, 'BR_stats5.mat'),'BR_stats5');
save(fullfile(path, 'BR_stats6.mat'),'BR_stats6');

% Write summary of image stats to summary .xlsx
countHeaders = ["Filename", "AvgCount_IF_mm2", "Inner_IF_mm2", "Outer_IF_mm2", "C1_IF", "C2_IF", "C3_IF", "C4_IF", "C5_IF", "C6_IF", "Placeholder",...
    "AvgCount_RM_mm2", "Inner_RM_mm2", "Outer_RM_mm2", "C1_RM", "C2_RM", "C3_RM", "C4_RM", "C5_RM", "C6_RM", "Placeholder2",];
areaHeaders = ["Filename", "AvgArea_IF_norm", "Inner_IF_norm", "Outer_IF_norm", "A1_IF", "A2_IF", "A3_IF", "A4_IF", "A5_IF", "A6_IF", "Placeholder",...
    "AvgArea_RM_norm", "Inner_RM_norm", "Outer_RM_norm", "A1_RM", "A2_RM", "A3_RM", "A4_RM", "A5_RM", "A6_RM", "Placeholder2",];
nnCountHeaders = ["Filename", "AvgnnCount_IF", "Inner_IF", "Outer_IF", "PercentAlone_IF", "nnC1_IF", "nnC2_IF", "nnC3_IF", "nnC4_IF", "nnC5_IF", "nnC6_IF",...
    "AvgnnCount_RM", "Inner_RM", "Outer_RM", "PercentAlone_RM", "nnC1_RM", "nnC2_RM", "nnC3_RM", "nnC4_RM", "nnC5_RM", "nnC6_RM"];
nnNearestHeaders = ["Filename", "AvgnearestDist_IF_mm", "Inner_IF_mm", "Outer_IF_mm", "nnD1_IF", "nnD2_IF", "nnD3_IF", "nnD4_IF", "nnD5_IF", "nnD6_IF", "Placeholder",...
    "AvgnearestDist_RM_mm", "Inner_RM_mm", "Outer_RM_mm", "nnD1_RM", "nnD2_RM", "nnD3_RM", "nnD4_RM", "nnD5_RM", "nnD6_RM", "Placeholder2"];
vimAreaHeaders = ["Filename", "AvgvArea_IF_norm", "vInner_IF_norm", "vOuter_IF_norm", "vA1_IF", "vA2_IF", "vA3_IF", "vA4_IF", "vA5_IF", "vA6_IF", "NSRatio_IF",...
    "AvgvArea_RM_norm", "vInner_RM_norm", "vOuter_RM_norm", "vA1_RM", "vA2_RM", "vA3_RM", "vA4_RM", "vA5_RM", "vA6_RM", "NSRatio_RM",];
areaSingleHeaders = ["Filename", "AvgASing_IF_norm", "ASingInner_IF_norm", "ASingOuter_IF_norm", "ASing1_IF", "ASing2_IF", "ASing3_IF", "ASing4_IF", "ASing5_IF", "ASing6_IF", "NSInner_IF",...
    "AvgASing_RM_norm", "ASingInner_RM_norm", "ASingOuter_RM_norm", "ASing1_RM", "ASing2_RM", "ASing3_RM", "ASing4_RM", "ASing5_RM", "ASing6_RM", "NSInner_RM",];
areaNetHeaders = ["Filename", "AvgANet_IF_norm", "ANetInner_IF_norm", "ANetOuter_IF_norm", "ANet1_IF", "ANet2_IF", "ANet3_IF", "ANet4_IF", "ANet5_IF", "ANet6_IF", "NSOuter_IF",...
    "AvgANet_RM_norm", "ANetInner_RM_norm", "ANetOuter_RM_norm", "ANet1_RM", "ANet2_RM", "ANet3_RM", "ANet4_RM", "ANet5_RM", "ANet6_RM", "NSOuter_RM",];

headers =  countHeaders;

% Create summary file if it doesn't exist
if isfile('ARA_QuantSummary.xlsx')
else
    writematrix(countHeaders,'ARA_QuantSummary.xlsx', 'Sheet', 'CountSummary');
    writematrix(areaHeaders,'ARA_QuantSummary.xlsx', 'Sheet', 'AreaSummary');
    writematrix(nnCountHeaders,'ARA_QuantSummary.xlsx', 'Sheet', 'nnCountSummary');
    writematrix(nnNearestHeaders,'ARA_QuantSummary.xlsx', 'Sheet', 'nnNearestSummary');
    writematrix(vimAreaHeaders,'ARA_QuantSummary.xlsx', 'Sheet', 'vimAreaSummary');
    writematrix(areaSingleHeaders,'ARA_QuantSummary.xlsx', 'Sheet', 'ASingleSummary');
    writematrix(areaNetHeaders,'ARA_QuantSummary.xlsx', 'Sheet', 'ANetSummary');
end

% Iterate through sheets and update with new summary martices
sheets = sheetnames('ARA_QuantSummary.xlsx');
sumLabel = array2table(cellstr(filename),'VariableNames', headers(1));

for i = 1:7
    if i == 2
        headers = areaHeaders;
    elseif i == 3
        headers = nnCountHeaders;
    elseif i == 4
        headers = nnNearestHeaders;
    elseif i == 5
        headers = vimAreaHeaders;
    elseif i == 6
        headers = areaSingleHeaders;
    elseif i == 7
        headers = areaNetHeaders;
    end
    mat = readtable('ARA_QuantSummary.xlsx', 'Sheet', sheets{i}, 'ReadVariableNames', 1); % Read in existing lines
    sumTable = [sumLabel, array2table([BIF_sum(i,:), BR_sum(i,:)],'VariableNames', headers(2:end))];
    matSum = [mat;sumTable]; % Merge tables together
    writetable(matSum,'ARA_QuantSummary.xlsx', 'Sheet', sheets{i});
end

disp('Data exported to ARA_QuantSummary.xlsx');
% End of ARA_Analyze.m